package math;

import java.util.Random;

public class NumberTheory {

	static final Random r = new Random();
	public long x, y;
	public int prime[], check[], size, upper = 32000;

	public void prime_table() {
		int i, j;
		size = 0;
		prime = new int[3432];
		check = new int[upper / 32 + 1];
		for (i = 2; i < upper; i++) {
			if (!bitAt(check[i / 32], i % 32)) {
				prime[size++] = i;
				j = i * 2;
				while (j < upper) {
					check[j / 32] |= 1 << (j % 32);
					j += i;
				}
			}
		}
	}

	/**
	 * @param n
	 * @return a factor of n
	 */
	public long pollard_rho(long n) {
		long d, x, y, rho;
		x = random(n);
		y = x;
		rho = random(n - 1) + 1;
		int i = 1, k = 2, max = (int) Math.pow(n, 0.25);
		max = max > 100 ? max : 100;
		while (i++ < max) {
			x = multiply_mod(x, x, n);
			x = (x + rho) % n;
			d = gcd(n + y - x, n);
			if (d != 1 && d != n) {
				return d;
			}
			if (i == k) {
				y = x;
				k += k;
			}
		}
		return pollard_rho(n);
	}

	public long gcd(long a, long b) {
		long c;
		while (b != 0) {
			c = b;
			b = a % b;
			a = c;
		}
		return a;
	}

	/**
	 * @param a
	 * @param b
	 * @return gcd(a,b), a*x+b*y=gcd(a,b)
	 */
	public long extended_euclid(long a, long b) {
		if (b == 0) {
			x = 1;
			y = 0;
			return a;
		}
		long t, d;
		d = extended_euclid(b, a % b);
		t = x;
		x = y;
		y = t - (a / b) * y;
		return d;
	}

	/**
	 * @param n
	 *            n is an odd integer greater than 2
	 * @param s
	 * @return whether n is a prime number
	 */
	public boolean miller_rabin(long n, int s) {
		long a;
		while (s-- > 0) {
			a = random(n - 1) + 1;// 1<=a<=n-1
			if (witness(a, n))
				return false;
		}
		return true;
	}

	/**
	 * @param a
	 * @param n
	 *            n-1=2^t*u
	 * @return
	 */
	private boolean witness(long a, long n) {
		int t = 0;
		long u = n - 1, x, tmp;
		while (u % 2 == 0) {
			u >>= 1;
			t++;
		}
		x = modular_exp(a, u, n);
		while (t-- > 0) {
			tmp = x;
			x = multiply_mod(x, x, n);
			if (x == 1 && tmp != 1 && tmp != n - 1) {
				return true;
			}
		}
		if (x != 1) {
			return true;
		}
		return false;
	}

	/**
	 * ax=b(mod n)
	 */
	public long modular_linear_equation_solver(long a, long b, long n) {
		long ans, i, d;
		d = extended_euclid(a, n);
		if (b % d == 0) {
			ans = x * (b / d);
			ans %= n / d;
			ans += n / d;
			ans %= n / d;
			for (i = 0; i < d; i++)
				System.out.println("The " + (i + 1) + "th answer is : "
						+ (ans + i * (n / d)));
			return ans;
		} else {
			System.out.println("no solutions");
			return -1;
		}
	}

	/**
	 * x=ai(mod mi) (0<=i<n)
	 */
	public long modular_equations_solver(int n, long a[], long m[]) {
		long d, lcm;
		int i;
		for (i = 1; i < n; i++) {
			d = extended_euclid(m[i], m[i - 1]);
			lcm = m[i] / d * m[i - 1];
			if ((a[i] - a[i - 1]) % d != 0) {
				System.out.println("no solutions");
				return -1;
			}
			extended_euclid(m[i] / d, m[i - 1] / d);
			a[i] += (m[i] / d * x * (a[i - 1] - a[i])) % lcm;
			a[i] = (a[i] % lcm + lcm) % lcm;
			m[i] = lcm;
		}
		return a[n - 1];
	}

	/**
	 * @param a
	 * @param b
	 * @param n
	 * @return a^b%n
	 */
	public long modular_exp(long a, long b, long n) {
		long d = 1;
		int len = bitLength(b);
		while (len-- > 0) {
			d = multiply_mod(d, d, n);
			if (bitAt(b, len)) {
				d = multiply_mod(d, a, n);
			}
		}
		return d;
	}

	/**
	 * @param x
	 * @param y
	 * @param n
	 * @return (x*y)%n
	 * 
	 *         x,y>=0
	 */
	public long multiply_mod(long x, long y, long n) {
		if (y == 0)
			return 0;
		return ((multiply_mod(x, y >> 1, n) * 2) % n + (x * (y % 2)) % n) % n;
	}

	public int bitLength(long n) {
		int l = 0;
		while (n > 0) {
			l++;
			n >>= 1;
		}
		return l;
	}

	public boolean bitAt(long n, int k) {
		return (n & (1l << k)) != 0;
	}

	/**
	 * @param max
	 * @return [0,max)
	 */
	public long random(long max) {
		return Math.abs(r.nextLong()) % max;
	}

}
